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Abstract. An inverse method is developed to determine the star formation history, the age-metallicity relation, 
and the IMF slope from a colour-magnitude diagram. The method is applied to the Hipparcos HR diagram. We 
found that the thin disk of our Galaxy shows a peak of stellar formation 1.6 Gyr ago. The stars close to the Sun 
have a solar metallicity and a mean IMF index equal to 3.2. However, the model and the evolutionary tracks do 
not correctly reproduce the horizontal giant branch. 
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1. Introduction 

The Hertzsprung-Russell diagram or its directly- 
observable counterpart, the colour-magnitude diagram 
(CMD) contains information about the age distribution 
of a stellar population, the initial mass function (IMF), 
and the metallicity of the stars. The distribution of the 
stars in such a diagram is essential for our understanding 
of the formation and evolution of the galaxy to which the 
stars belong. 

In order to infer the star formation history (SFH), 
one usually tries to match the CMD as closely as pos- 
sible by model populations computed with various sce- 
narios for the SFH. The quality of the agreement between 
observed and predicted CMDs is assessed either qualita- 
tively by visual inspection or quantitatively, e.g. by the 
X 2 -test (cf. Dolphin, 1997). Often, one represents the SFH 
by a number of discrete values or by a suitable analytic 
form, and adopts either a fixed metallicity or a certain 
age-metallicity relation. This approach, which could be 
called the direct or synthetic method, has been widely 
used (Tolstoy et al. 1993, Gallart et al. 1996, Haywood et 
al. 1997). Since one usually keeps the number of free pa- 
rameters for the underlying model as small as necessary 
or practical, the adopted functional forms for the SFH 
and age-metallicity relation (AMR) impose certain limi- 
tations on what type of model populations are possible 
to be considered. Thus the obtained best fit is only opti- 
mal in the context of the adopted model, which may not 
necessarily represent the real galaxy. In what sense the 



interpretation will be limited by these constraints is diffi- 
cult to assess accurately, and can only be done by experi- 
mentation. In most recent work (Harris & Zaritsky 2001, 
Dolphin 2001) these limitations are being overcome by us- 
ing more efficient minimization techniques. The technique 
of inverse methods (Craig & Brown, 1986; Tarantola & 
Valette, 1982a, b; Twomey, 1977) deals with such a prob- 
lem in the opposite way: One tries to determine the func- 
tional form of e.g. the SFH with as much freedom as pos- 
sible - with a resolution that depends on the informa- 
tion contained in the data - under the constraint that it 
matches the observed data. Such a method has been ap- 
plied to the CMD by Hernandez et al. (1999): this work 
presents a non-parametric method for the maximum like- 
lihood solution of the SFH, through the iterative solution 
of an integro-differential equation. The method presented 
here is similar and has the advantage of determining not 
only the SFH from CMD but also the IMF slope and 
the AMR with a temporal resolution as high as possi- 
ble. Specific tools are presented in order to estimate the 
validity of the inverse procedure like the a posteriori co- 
variance and the resolution. We apply our method to the 
Hipparcos data for solar neighbourhood stars in order to 
give constraints on the Galaxy thin disk formation. 



2. Relation between SFH and the HR diagram 

The distribution of stars in a CMD is determined by the 
following ingredients: 
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— the SFH rb(t) specifies the number of stars of all masses 
born as a function of age t; 
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— the initial mass function (IMF) gives the number N 
of stars in each generation per unit interval in stellar 
mass M. A convenient form is a power-law 

dN a: M~ T dM 

where the index T = 2.35 refers to Salpeter's (1955) 
function. As usual, it is assumed to be independent of 
time; 

— due to the chemical evolution of a galaxy, the chemi- 
cal composition of the gas from which stars are born 
changes with time. This is described by an age-metal- 
icity relation (AMR) Z(t); 

— the results of stellar evolution calculations give the 
properties - temperature, luminosity - of a star of 
given mass and metallicity for any given time after 
its birth; 

— stellar atmosphere models give the photometric colours 
for any star. 

For comparison with models, one bins the CMD in 
some suitable way in colour index and magnitude. Then 
the number density of stars in a given bin q can be written 
as 

D q = J iP(t)F q (t,Z(t),T)dt (1) 

the superposition of the contributions from all stellar pop- 
ulations of single age t, metallicity Z, and IMF slope T. 
Since the densities F q (t, Z, T) would be computed from 
stellar evolutionary tracks and the photometric colours, 
the interpretation of the observations D q — D q (obs) con- 
sists in solving Eqn. [j] for the SFH tp. 

In principle, one would like to determine all three func- 
tions SFH, IMF, and AMR, i.e. functions defined on an 
infinite number of points. However huge the observed set 
of data might be, it remains finite. Thus, even if there 
were no noise and observational uncertainties, one can de- 
rive the three functions only up to a certain limit of detail. 

The usual approach is to assume that the functions 
have some specified analytical form described by a few 
parameters - such as the exponent in the power-law for the 
IMF, or a simple SFH with a few rectangular starbursts - 
and the objective is simply to find the best model of this 
type that reproduces the observations. In doing so, one 
does not know of what (time) resolution could be reached 
with the data available, and thus one may not make full 
use of what the observations could furnish. 

3. The inverse method 

The determination of the full functional form of the SFH 
(or AMR or IMF) from a CMD is an under-determined 
problem, because the amount of observational data is only 
finite and thus cannot provide the information for every 
detail of the function. Application of a straight inversion 
technique to Eqn. ^ could be very sensitive to the noise in 
the data, and could well give mathematically correct but 
unphysical results (Craig & Brown 1986). The problem 



can be regularized by demanding that the true solution is 
smooth in some suitably quantified sense (Twomey 1977, 
Craig & Brown 1986). Thus one reduces the number of 
free parameters, and the problem becomes well-posed. 

Tarantola & Valette (1982a) use a Bayesian approach 
which describes how the a priori knowledge about the 
functions - which might be the null information - is 
changed by the information contained in the observational 
data. If the data contain no information about a param- 
eter, one would merely retrieve its a priori value. The a 
posteriori probability density f pos %(M\D) for the vector 
M containing the unknown model parameters, given the 
observed data D, is linked by Bayes' theorem 

f post (M\D) <x L(D\M) ■ f prioI (M) 

to the likelihood function L and the prior density function 
for the parameter vector. The factor of proportionality is 
obtained by normalization J f post (M\D)dM = 1 over all 
parameter space. 

The theoretical model shall be described by an oper- 
ator g which connects the model parameters M with the 
predicted data 

-^predicted = 9{M) 

which is to agree as closely as possible with the observed 
data. If we assume that both the prior probability and the 
errors in the data are distributed as Gaussian functions, 
the posterior distribution becomes 

/ post (M|£>) oc eM~\{D - g(M)) T C^(D - g{M)) 

-l(M-M(0)) T C^(M-M(0))) 

The prior is specified by the value M(0) of the parameter 
vector and its variance-covariance matrix Co- The matrix 
Cd specifies how the variances in the observed data are 
obtained (e.g. if a non-Gaussian distribution of the obser- 
vational errors is to be approximated) and whether one 
needs to take into account correlations between the indi- 
vidual data. 

In this method the model parameters may include sin- 
gle value parameters as well as entire functions. For simple 
parameters, Cq is the matrix of its a priori variances and 
a priori covariances between each parameter. The vari- 
ances and covariances would be chosen to be large, if we 
do not have any initial knowledge about the parameter 
values. On the other hand, if the values are accurately 
known, one would use variances correspondingly smaller. 
For a function, Co is a functional operator which has for 
its kernel the auto-correlation function. Most commonly 
one uses a Gaussian kernel: 

C (x,x') = <T (x)<T (x')exp y- — ^ - j 

ctq(x) - which can be any function of x — describes 
how strongly the solution is allowed to fluctuate between 
points, and £o is the smoothing length for the solution. 
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Other types of kernels may be employed, e.g. an exponen- 
tial kernel which would smooth the solution differently. 

The best estimator M for M is the most probable value 
of M, knowing the set of data D. This condition is reached 
by minimizing the quantity 



-{D-g(M)) T C»\D-g{M)) 



+ -(M - M(0)y C^iM - M(0)) 



(2) 



which corresponds to a maximum likelihood condition. 

Since the operator g is non-linear, the solution for the 
parameter vector must be done iteratively (Tarantola & 
Valette 1982a, 1982b): 

M(k + l) = M(0) + C o G T (k)(C D + G(k)C o G T (k)y 1 

(D + G(k)[M(k) - M(0)] - g(M(k))) (3) 

with k counting the number of iterations, and the matrix 
of partial derivatives 



G(k) 



dg(M(k)) 
DM 



For use in subsequent equations we shall abbreviate: 
S :=C D + G(k)-C Q -G T (k) 

3.1. Formulation of the problem 

The CMD is binned both in M v and (B-V), by divid- 
ing the range of interest of My-values into m intervals 
of equal length AMy. Likewise, the range of (B-V) val- 
ues is divided into n intervals of equal length A (B-V). 
Then, D(i,j) is the number of stars having their magni- 
tude My between My 4 and My,j + AMy as well as their 
colour index (B-V) between (B-V), and (B-V),- + A(B- 
V). For a less cumbersome notation, we use a single index 
q = (i — l)n + j = 1 . . . nm corresponding to bin 

Synthetically-generated stellar populations for single 
ages t, metallicities Z, and IMF slopes T are used to spec- 
ify which fraction of the total number of stars is found in 
each bin q 



F q (t, Z, T) AMy A(B-V) 



(4) 



Convolving these populations with the SFH by Eqn. [I] 
gives the number density D s of stars in that bin. 

Since the star formation rate is always positive, chang- 
ing over from ip to a(t) 



a(t) := lnty>(t)M,) 



(5) 



where ipo is a constant, will enforce this physical fact. We 
shall assume that the uncertainties in the function a(t) 
follow a Gaussian law. This implies that the errors in ip(t) 
follow a log-normal distribution. 

Since the position of isochrones in the CMD is a nearly 
logarithmic function of the age, it is more appropriate 
to use u := lg(i) instead of the linear age t. With dt = 



10" La(10)du and the abbreviation h(u) := 10" ln(10), 
Eqn. is written as 



D„ 



ipo exp(a(u))F q (u, Z[u), T)h(u)du 



g q (M) 



(6) 



This defines g q as a non-linear operator acting on the pa- 
rameter vector M which consists of the unknown functions 
a{t) and Z(t), and the unknown parameter T: 

We note in passing that the IMF could of course also be 
taken as an unknown function cf>(m). This would merely 
increase greatly the dimension of parameter space, result- 
ing for the numerical solution of the equations in greater 
memory requirements and unacceptable execution times. 
Furthermore, an arbitrary IMF would be largely degener- 
ate with the SFR, at least on the main sequence. 

3.2. The base models 

The operator g q (Eqn. ^) generates the model CMDs from 
a set of base models F q (u,Z,T) (Eqn. 0) by integrating 
over all ages. The base models are computed from syn- 
thetic stellar populations of a single age t, metallicity Z, 
and specified IMF; they give the probability of finding a 
star in bin q in the CMD, depending on the model param- 
eters, but we also take into account observational selection 
effects: 

— synthetic theoretical stellar populations are computed 
from the isochrones of Bertelli et al. (1994) which di- 
rectly give the positions in the (My, (B-V))-diagram 
for stars in the mass range 0.6 to 120 M Q ; 

— the IMF is a power-law between the stellar mass limits 
of 0.6 and 80 M Q ; 

— the stars are uniformly and randomly distributed in 
space; 

— the limiting magnitude is 8.0 mag in V band. 

A grid of models was computed, each comprising 2 mil- 
lion stars (before applying selections), for any combination 
among 50 ages distributed linearly between u = 6.6 and 
10.3, 10 metallicities between Z = 0.0004 and 0.05, and 5 
IMF slopes between 2.0 and 4.5. This was done through 
linear interpolation in age and metallicity from the Padua 
isochrones. 

All integrals are evaluated numerically, with the rect- 
angle method being sufficiently accurate. 



4. Determination of SFH and IMF slope 

We shall first consider the case of deducing simultaneously 
the SFH and the slope of the IMF: The parameter vector 
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M consists of the unknown function a(t) and the unknown 
parameter T: 

' a(*) 

r 



M = 



The metallicity is considered to be known and constant. 
Then, the matrix of partial derivatives is 



/ dgi_ dgi_ \ 

I da ar \ 



G(k) 



da ar 



\ da ar 



with s = nm the total number of data bins. 

Equation ^| gives the procedure to improve an esti- 
mate of the parameters. For better legibility, we drop the 
dependence on the iteration number k of G, a(u), and T. 
The (z,j)-th component of the matrix product is (for all 
i,j = l,...s) 



(GCoG T )ij 



C a (u, u')i() exp(a(u) + a(u')) 
x Fi (u, T)Fj (vf, T)h{u)h{u')dudu' 
ipo exp(a(«)) h(u)du 



tp exp(a(u)) 



dT 

dFj(u, r) 

dr 



h(u)du 



The variance-covariance matrix Co of the prior distribu- 
tion of the parameters is 



Cn 



C Q 
al 



and the covariance function C a 

n I l\ 2 f ( U ~ U ') 2 

C a (u, u ) = a a exp 



(7) 



with £ a the correlation length in lg(t). 

The i-th component of the vector V(k) = D — g(M) + 
G(k)(M — M (0)), which depends on the iteration index k, 
is computed from 



V. ^ H + / Voexp(a(u))(a(w) - l)F t (u,T)h(u)du 
Tp exp(a(u))(r - T ) dFl ^ r ^ h(u)du 

Defining the vector W(k) = S~ 1 V(k), the fc+l-st estimate 
for the parameters is computed from their fc-th values by 

a k+1 (u) = Wj(k) J C a (u,u')^ exp(a k {u')) 

i 

xFi(u', Tk)h(u )du 

and 

r fc+1 = T + Wi(k)o$ f ^ exp(a fc ( U ')) 

i 

dF(u\T k ) 



dT 



-h(u')du' 



The estimator for ip from a is given by the relation (e.g. 
Saporta, 1990): 



ijj(u) = V>o exp a(u) + 



'a(u) 



with the dispersion of the posterior distribution for a(t) 

o>(„) = V'oy / 'exp(2Q ! (u) + a 2 a(u) ) (exp(al (u) ) - l) 

where cr a ( u ) is the dispersion of the posterior distribution 
of a (Eqn. |9| below) . 

In the computations, the derivatives such as the matrix 
G and dF/dT are evaluated numerically. The functions 
OF/dT are interpolated by a 4 th order polynomial. 

4.1. The posterior variance-covariance matrix 

The internal errors made by the method on the estimated 
parameter values can approximately be computed. This is 
done by a second-order expansion of the posterior density 
distribution in the neighbourhood of the best solution: 

G £,[ — Cq — CqG t S """GCo (8) 

For the simple parameter of the IMF slope one gets 



<r f = <t t ^/(1-G t S- 1 G), 
and for the function a(u) : 



<?a(u) = A/ ~ C a G T S 1 GC Q 



4.2. Resolving kernel and mean index 



(9) 



Observational errors, like the finite photometric accuracy, 
will degrade the information contained in the CMD, and 
thus cause an additional loss of resolution in age. The 
concept of the resolving kernel, introduced by Backus & 
Gilbert (1970), permits us to compute how much a stellar 
population for a given age is spread out in age. Suppose 
that we knew the true model M true , then the observed 
data would be 

D = g(M tTue ). 

If the operator g were linear, so that G k = G = g, then 
Eqn. |^ shows how the deviation of the true model param- 
eters from the initial guess is expected to be degraded into 
the observed one 

M - Mo = C G*S- 1 G(M tTUC - M ). (10) 

The operator - called the resolving kernel - 

K(u,u') :=C G*S- 1 G 

describes this degradation of information. Though it ap- 
plies strictly for linear models only, it remains a good ap- 
proximation in our problem as well: the density of stars in 
a bin is a linear function of the SFH. Because the values 
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for the SFH remain within a factor of 10 or so, the loga- 
rithmic relation with the parameter a does not represent 
a strong non-linearity. 

The relation between the resolving kernel K and the 
a posteriori variance-covariance operator is computed 
as follows from Eqn. ||: 



C 



M 



(I - C G T S' 1 G)C = (J — K)C 



this yields: 



K — I C^Cq 1 



I is the Dirac function if the unknown parameter is a func- 
tion (the SFH in our case). So, if the posterior variance- 
covariance operator increases, the amplitude of the resolv- 
ing kernel decreases and the width increases. 

Another important and useful concept is a measure of 
the information present in the data. This is closely linked 
to the resolving kernel: Suppose that the mean value of 
the SFH is within the width of the kernel K, so Eqn. [To| 
gives 

(M - M(Q))(«) = (M truc - M(0)) meail 
The integral is called the mean index I(u) 
I(u) := / K{u : u')du' 



K(u,u')du'. 



which has the following meaning: If I(u) has very low val- 
ues, the expected model would be close to the initial guess, 
irrespective of how much the initial model would differ 
from the true model. The quality of the SFH estimator 
thus is poor. But if I(u) « 1 the deduced model would be 
close to the true model. 

We note that if we take an a priori variance-covariance 
operator with too large a smoothing length £ Q , we will 
obtain a good mean index. However, the a posteriori res- 
olution will be poor. 

5. Simultaneous derivation of SFH and AMR 

We now address the problem of determining simultane- 
ously two functions, the SFH and the AMR. The param- 
eter vector is 

a(t) 
Z(t) 

The variance-covariance matrix of the parameter priors is 
taken as: 

C a 
C z 



M 



Co 



i.e. we assume that SFH and AMR are a priori uncorre- 
cted. This is done for simplicity. We note that chemical 
evolution of the galaxy would impose a strong relation 
between the two functions and could be used to compute 
such a correlation. 

C Q is given by Eqn. [j], and Cz is given by 



C z (u,u') = erf exp 



(u- 



u') 2 



f 2 



The derivation is nearly identical as before, with a(u) and 
Z(u) depending on the iteration index k: 



(GCoG*)i tj = J C a (u,u')<eMa(u) + a(u')) 

xFi(u, Z{u))Fj{u\ Z(u'))h{u)h{u')dudu' 
+ [ [ C z (u,u')^ 2 exp(Z(u) + Z(u')) 



dZ 



dZ 



In addition to the iteration for the SFH which proceeds 
as before, we have for the AMR 

Z k+1 {u) = Z Q {u) + Y,W l {k) J C z (u,u')i) exp{Zk{u')) 

i 

xFi(u',Z k (u'))h(u')du' 
where Zo(t) is the a priori age-metallicity-relation. 

6. Validation of the method by simulations. 

In practice, the large amount of computing time necessary 
makes it prohibitive to determine simultaneously the SFH, 
AMR, and the IMF slope. Therefore, in the following tests 
we show an inversion of SFH and IMF slope only. Also, to 
limit the number of bins and hence the computing time, we 
consider only the histogram of the (B-V)- colour index. 
The colour of the turn-off is very sensitive to the age of a 
stellar population, thus one obtains a good resolution of 
the SFH. The range from —0.3 to 1.7 is divided into bins 
of 0.02 mag width, the observational error. It is not useful 
to use cell sizes smaller than this error. 

The variance-covariance matrix C d for the data is as- 
sumed to be diagonal, as we shall assume that there are no 
correlation between the numbers of stars in different bins. 
To approximate the Poisson noise in the stellar numbers, 
we take the variance equal to the number of stars in a bin. 

Since our method is an iterative procedure, it yields 
models that approach the data successively closer with 
each step. The distance between model prediction and ob- 
served data can be measured by the value of \ 2 : 



E 



{g s (M)-D s )< 



The convergence is usually monotonic. Typically 10 iter- 
ations suffice to achieve a value of x 2 stable within 10 -2 
for the determination of the SFH. 

In our formalism, the errors in the data are assumed 
to be normally distributed. This is an approximation valid 
only if each bin contains a large number of stars, whereas 
in reality the errors in the histogram are distributed with 
a Poisson law. 

As simulated observational data we took a stellar pop- 
ulation model having a fixed (solar) metallicity, a power- 
law IMF with slope 2.25, and a star formation history in- 
volving several bursts at different ages and with different 
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strengths (shown as the full line in the middle left panel 
of Fig. |). 

We have to specify the parameter a a which quanti- 
fies how large a fluctuation we accept in the SFH. We 
chose it in such a way that one obtains reduced \ 2 = 
0.5. ..1.5. Smaller values would show that the data are over- 
interpreted, but larger values indicate that the model re- 
produces the observations only poorly. It turns out that 
a a = 1 is a good choice. The value of the smoothing length 
£ Q is normalised as er^ ~ cste/£ Q . 

In Fig. [l] we show the results obtained for the data 
to which Gaussian noise had been added with an ampli- 
tude of ctmv = 0.3 ma g an d cr (B-v) = 0-01 mag. The best 
solution has a reduced x 2 = 0.77, and as the top pan- 
els show, the stellar population is closely reproduced. The 
slope of the IMF is T = 2.16 ± 0.08, i.e. slightly flatter 
than had been assumed. The middle panels depict that 
the deduced SFH well matches all the features of the as- 
sumed SFH older than about 100 Myrs. One notes that 
even the small SFH feature near 1 Gyr is reproduced. 
However, the young populations are not very well rep- 
resented because the number of corresponding 'observed' 
stars is rather small (< 100) and, in addition, the B — V 
colour of the main sequence becomes insensitive to age 
for young ages. This is also evident in the bottom panels: 
while the resolving kernel shows that the relative resolu- 
tion in age is rather good and nearly the same for ages 
larger than 100 Myrs, the resolution (in lg t) deteriorates 
for younger populations. Finally, the level of information 
in the CMD - as measured by the mean index - drops 
below about 10 Myrs. 

Figure || presents the results for the same 'observa- 
tions', but degraded with larger values <tm v — 0.5 mag 
and CT(b-v) = 0.05 mag, which is evident from the CMDs. 
The reduced x 2 obtained here is 0.64. As before, the IMF 
comes out somewhat too flat with V = 1.99±0.12. The in- 
creased noise level in the data destroys all the information 
about the details of the SFH for ages up to about 1 Gyr; 
only the oldest starbursts can be detected. The time res- 
olution of the method has become much worse than Fig. 
|l| for all ages except around 1 Gyr; the drop of the mean 
index now occurs at about 30 Myrs. 

Figure || presents the results of the SFH and AMR in- 
version from a simulated CMD composed of an old metal- 
poor population (Z = 0.004) and a young metal-rich one 
(Z = 0.02). Because the zero age main sequences are 
shifted with metallicity, the SFH- AMR information is not 
degenerated in the colour-magnitude diagram. The prior 
taken for the inverse process is tpo = 1, a(t) = and 
Z{t) = 0.01, implying constant SFH and AMR. The con- 
vergence occurs after about 20 iterations. 

These test results demonstrate that the method is ca- 
pable of deducing reliably the history of star formation 
in rather fine detail. It also permits one to assess quanti- 
tatively the time resolution that is possible with a given 
observational data set, as well as the information it con- 
tains. 



7. The local CMD from Hipparcos data 

Hipparcos gives astrometry and photometry for stars 
within about 1 kpc of the Sun. Most of these stars belong 
to the thin disk, only about 5 % of them being members of 
the thick disk and halo (Robin et al., 1996). This sample 
gives us valuable information about the evolution of the 
thin disk of our Galaxy. This is not limited to the solar 
neighbourhood proper, but pertains also to a larger por- 
tion of the disk: because of the radial diffusion of stars in 
the disk (Wielen, 1977), one finds presently, near to the 
Sun, stars of different metallicities and ages. 

7.1. The Hipparcos data 

We have selected a sample of 13520 stars from the 
Hipparcos catalogue (Perryman et al.; ESA 1997) which 
meet these criteria: 

— single stars 

— the apparent magnitude limit V depends on the lati- 
tude b : 

V < 7.3 + l.l|sin(6)| 

for all spectral types; 

— observed parallaxes should be larger than 5 mas. 

Double stars may cause some problem: confirmed dou- 
ble stars in the Hipparcos catalogue comprise about 20% 
of all stars. Also, it is well known that among a stellar 
population about 50 % are double stars. Hence, we under- 
estimate the presence of double stars, which could cause a 
bias in the following way: binaries will appear to be located 
on the red side of the main sequence. They will be inter- 
preted by the model as being cither of higher metallicity 
or in a more evolved state. Thus, the derived populations 
could be too old or too metal-rich. 

The choice of the bin sizes is determined by the obser- 
vational uncertainties of the data: For the error on the par- 
allaxes 7r we take a constant value as a reasonable approxi- 
mation a n = l(mas). This results in an uncertainty of the 
absolute magnitudes of ctm v — 2.17 a n /n which reaches 
0.5 mag at 200 pc. The error in the colour (B-V) is taken 
as constant at 0.02 (mean error in the Hipparcos catalog). 
Therefore, we divide the HR diagram into cells of 0.5 mag 
in My and 0.02 in (B-V). To prevent nearly empty cells 
with their large relative fluctuations of occupation num- 
bers to influence the results too strongly, we increase the 
cell size, if necessary, so that each bin shall contain at least 
5 stars. The resulting cell layout is presented in Fig. |[ 

With the limit on apparent magnitude and the local 
density of the stars, one expects only very few stars below 
0.6 M . Thus, the isochrones of Bertelli et al. (1994) which 
pertain only to stars above 0.6 Mq are entirely adequate. 

Absolute magnitudes and colour indices are corrected 
for interstellar extinction with the model of Vergely (1998) 
which gives the opacity at each point of the space in the 
250 pc sphere surrounding the Sun, derived with a to- 
mographic technique from Hipparcos data and Stromgren 
photometry. 
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Fig. 1. The top panels show the simulated 'observed' CMD (left) with ctm v = 0-3 and C(b-v) = 0.01 and the best 
model (right). The middle left panel shows the assumed SFH (dashes) and the best model (full line), on the right is the 
best model at ±lu. The bottom left panel shows the resolving kernel K(t,t') as a function of age t for log(t') = 7.3, 
8.1, 8.8, and 9.6 (full, dashes, dot-dashes, and dots), and at right the mean index. Note that the fit for the tests is 
performed on the B — V distribution only, not on the full 2-D CMD distribution. 



7.2. The base models F(u) 

The construction of the base models follows that described 
in Sect. 3.2, but we now take into account the above- 



mentioned selection criteria for the observational data. 

Furthermore, for solar neighbourhood stars, we must 
also take into account that older stars are distributed up 
to greater heights above the Galactic plane than younger 
ones, because the amplitude of their oscillatory motions 
perpendicular to the plane increases with age (cf. Wielen 
1977). 

In this study, we restrict our SFH-AMR analysis to 
the thin disk stars. The contamination by thick disk pop- 
ulation stars can be neglected in a first approach, because 



only about 5% of stars belong to the thick disk in the 
Hipparcos sample. The age-dispersion relation given by 
Gomez et al. (1997) is valid only for thin disk population 
stars and shows clearly that there is no significant dynamic 
evolution of the thin disk, after 4-5 Gyr. More accurately, 
the age-dispersion relation is given for old thin disk stars 
with age < 9 Gyr and saturated at 15-17 km/s. 

The relation between age t ([Gyr]) and dispersion of 
the vertical velocity component W ([km/s]) has recently 
been derived from Hipparcos data (Gomez et al., 1997): 

o w {t) = 17.0 - 12.0cxp(-</2) 

Assuming isothermality, the density of coeval stars of age 
t is constrained by the vertical potential <j)(z) at height z 
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Fig. 2. Same as Fig. || but with ctm v = 0.5 and d( B -v) = 0-05 



as: 



p(t, z) = p(t, 0) exp 



,{t) 



where p(t, 0) is the density at z = and 4>(z) is the po- 
tential. 

Recent determinations of the total local mass density 
p , based on Hipparcos data, give p from 0.076 M Q pc~ 3 
(Creze et al, 1997,1998) to 0.10 M Q p C - 3 (Holmberg & 
Flynn, 2000). The surface mass density is determined from 
measurements of the force at larger distances (0.5 to 1 kpc) 
from the galactic plane: recent determinations of £o range 
from 48 to 56 M pcr 2 , (Kuijken & Gilmore 1989, Flynn & 
Fuchs 1994). These two independent sets of constraints on 
the vertical potential may be combined in a single formula: 



$(z) = 0.027 E (\/z 2 +D 2 - D 



Pcft 



(11) 



Here, we will use the coefficients D — 240 pc, So = 
48M©pc -2 and p c g = 0.0105 M Q pc~ 3 , this implies a lo- 
cal mass density of 0.081 M Q pc~ 3 and a disc surface mass 
density of 48 M Q pc~ 2 . 

7.3. Determination of IMF slope and SFH 

We first consider the metallicity to be constant, at solar 
value Z — 0.02. The resultant fit is not good (reduced 
X 2 = 4.9), so that the model is to be rejected. 

The derived IMF slope is T — 3.2 ± 0.1, steeper than 
Salpeter's, but quite close to the slope of 2.7 found by e.g. 
Kroupa, Tout & Gilmore (1993) in the solar neighbour- 
hood for stars above 1 M Q . It is these stars that contribute 
to most of the data in our sample. 

Main sequence magnitudes lie between about 6 and 
—2 mag, corresponding to masses between 0.6 and 8 M Q , 
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Fig. 3. The top panels show the simulated 'observed' CMD (left) with au Y = 0.3 and crm-v) = 0.01 and the data 
fitting in number of stars per bin (right). The middle left panel shows the assumed SFH (dashes) and the best model 
(full line); on the right is the assumed AMR (dashes) and the best model (full line). The bottom left panel shows the 



resolving kernel K (t, t') as a function of age t for log(f') 
and at right the mean index 



7.3, 8.1, 8.8, and 9.6 (full, dashes, dot-dashes, and dots), 



and lifetimes of greater than 40 Myrs. Thus, information 
about the last 40 Myrs is not present. This is seen in the 
turn-down of the mean index (cf. Fig. §). The prior taken 
for the inverse process is (see Eqn. |5|) a(t) = 0, implying 
a constant SFH. 

The deduced SFH is presented in Fig. ^. The resolution 
is rather poor; from the widths of the kernels one gets 
0.4 dex in the relative age. The a posteriori error is quite 
large for the youngest populations and the results are not 
significant for ages below 30 Myrs, as clearly shown by the 
low mean index which drops at about 50 Myrs. 



7.4. Simultaneous determination of SFH and AMR 

As the model with constant metallicity has to be rejected, 
we now investigate models with a fixed IMF slope of 
r = 3.0, but leaving both the SFH and AMR as free pa- 
rameters. As priors for the inverse process we take a(t) = 
and Z(t) = 0.02, i.e. constant SFH and AMR. 

If we assign equal weights to main sequence and giant 
stars, we obtain a rather poor value of reduced \ 2 — 3.7. 
The model population does not fit the main sequence, in 
particular below it, many metal-poor stars are predicted. 
The reason is that the prominent red clump - composed 
of horizontal branch stars of nearly solar metallicity - con- 
tains many relatively blue stars which in the context of the 
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Fig. 4. The colour-magnitude diagram of the single stars 
of the Hipparcos catalogue, with the adopted binning in 
both dimensions 



presently-used stellar tracks correspond to low metallicity 
isochrones. 

If we reduce the weight for the giants by a factor of 
two, a rather satisfactory fit with reduced \ 2 — 1-9 is ob- 
tained. Figure || shows that the main sequence is well re- 
produced except for the odd bin, but that a major discrep- 
ancy remains with the red clump where the observed stars 
tend to be systematically bluer than the model. Very re- 
cently, Girardi et al. (2000) published updated isochrones. 
Detailed comparison with the isochrones of Bertelli et al. 
(1994) reveals that the horizontal branch is shifted to the 
blue by about 0.05 mag in (B-V) which is the right amount 
to remove the discrepancy (as shown by test calculations). 

However, this shift could be explained by a dispersion 
in the AMR (not taken into account in this study) or the 
presence of unresolved giant binary stars. 

7.4.1. The star formation history 

The SFHs obtained from the two approaches are quite 
similar, having a rather broad peak near an age of 1.6 Gyr. 
The SFH at the present time deduced with the constant 
metallicity model is rather low, while in the model with 
variable metallicity (Fig. [7]) it is almost as high as this 
peak. Note howeverthat the information available allows 
a reliable deduction of the SFH with good resolution only 
for ages larger than about 100 Myrs. The age resolution 
achieved in the second model is somewhat poorer, due to 



in 




500 1000 1500 

z (pc) 



Fig. 5. The vertical potential (continuous line) force given 
by Kuijken & Gilmore (1989) and the modified force 
(dashed line) used here in order to include the most re- 
cent determinations of a local mass density in the solar 
neighbourhood. 



the information now being used to derive the AMR as 
well. 

A first preliminary analysis of the Hipparcos CMD by 
Bertelli et al. (1997) yielded that star formation occurs 
between ages of 0.1 and 10 Gyrs, with indication of a dis- 
continuity at about 1.5 Gyr. Also, the broad red clump 
indicated a spread in metallicity. From a rather qualita- 
tive comparison they find that the SFH should have been 
constant or increasing from 10 to 1.5 Gyrs, and be re- 
duced by a factor 2 or 3 after that time. Our results fully 
confirm these findings, showing a broad and rather well- 
defined peak of star formation activity at about 1.5 Gyrs 
with the SFH being substantially lower before and after. 

Bertelli and Nasi (2000) use a comparable set of 
Hipparcos stars to analyse the SFH in the solar neighbor- 
hood. Their principal conclusion is that the SFH increases, 
in a broad sense, from the beginning to the present time; 
that is in agreement with the SFH obtained in this pa- 
per. They have tested different IMF slopes and preferred 
the IMF Salpeter slope in order to fit the main sequence. 
However they underline that an IMF slope of 3.35 is an 
acceptable solution, particulary to find the best ratio be- 
tween the number of stars in the He burning phase and 
that in the main sequence phase. Our analysis shows that 
this IMF slope can fit very well the main sequence, too. 
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Fig. 6. The star formation history of the solar neighbour- 
hood, assuming a constant stellar metallicity Z = Z©. 
From top to bottom: The derived SFH with standard er- 
rors, the resolving kernel K(t, t') as a function of age t for 
log(i') = 7.3, 8.1, 8.8, and 9.6 (full, dashes, dot-dashes, 
and dots) and the mean index. 



Hernandez et al. (2000) use an inverse method to de- 
rive the SFH over the last 3 Gyr. They show clearly a 
bump in the SFH at about 2 Gyr. Similar to their results 
we do not find the decreasing SFH activity between 1 and 
2 Gyr as found by Rocha-Pinto (1999) from chromospheric 
activity. 

The determinations of the IMF (cf. Scalo 1986) yield 
that the present SFR should be nearly equal to the past 
average SFR. Our study gives the same result: 



/ 10Gyr ^(t)dt/10Gyr 



0.9. 



A quite different approach can be obtained from the 
chemical evolution of the solar neighbourhood: Rocha- 
Pinto & Maciel (1997) derived from the metallicity distri- 
bution of the G-dwarfs and the AMR the SFH (cf. Eqn. |l3] 
below) taking into account the observational scatter. They 
demonstrate that for several observed AMRs one arrives 




og(age) 



Fig. 7. Similar to Fig]?], but for the simultaneous inversion 
for SFH and AMR. The second panel from the top shows 
the deduced AMR with standard errors 



at rather similar SFHs: a broad peak at about 8 Gyrs 
ago with about twice the present SFH and an epoch of 
low SFH (one half or less of the present value) between 
1 and 3 Gyrs ago. As this approach requires the slope of 
the already rather uncertain and noisy AMR, the results 
are quite sensitive to the adopted AMR, e.g. the mean 
AMR of Edvardsson et al. (1993) gives a nearly flat SFH. 
At an age of 1.5 Gyrs, when the Hipparcos CMD gives 
a peak SFH, Rocha-Pinto & Maciel rather find a mini- 
mum. Though their result seems statistically significant, 
one should emphasize that this approach can yield only 
a resolution which is nearly constant in (linear) time, be- 
cause the metallicity changes nearly linearly in time with 
some flattening near the present. On the contrary, the dis- 
tribution of stars in a CMD is sensitive to the logarithm 
of the age. Therefore we feel that the Hipparcos CMD 
gives a stronger constraint on the SFH in the past few 
Gyrs. Certainly, this discrepancy between chemical evolu- 
tion and stellar population needs further attention. 

As a further check we investigate whether our SFH is 
in agreement with the past gas consumption in the solar 
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Fig. 8. Comparison of model and observed CMD: The 
/ (\) indicate the bins whose number of model stars is 
smaller (greater) by at least 3a than the observations 



neighbourhood: in a closed-box model the evolution of the 
surface gas density S g is given by (Tinsley, 1980) : 

^L = -(l-R )m . (12) 

The constant (1 — R) is the lockcd-up mass fraction, which 
is about 0.8 for a Salpetcr IMF, and rather insensitive to 
the IMF. The total surface density of the mass present in 
the disk was derived from stellar dynamics by Kuijken & 
Gilmorc (1991) as 48 M0pc~ 2 . The current gas surface 
density is 6.6 Mope -2 estimated by Rana (1992) as the 
sum of molecular and atomic hydrogen. Wyse & Gilmore 
(1995) argue for a somewhat higher present day gas frac- 
tion of 25 percent. Using the total mass as the initial con- 
dition, we integrate Eqn. |lj. The evolution of the gas den- 
sity shown in Fig. || matches the present day gas fraction 
rather well. 

7.4.2. The age-metallicity relation 

In their preliminary study, Bertelli et al. (1997) empha- 
sized that the colour extension of the horizontal branch 
clump as well as the width of the main sequence pro- 
vide evidence for an extended range (Z = 0.008... 0.03) 
of metallicities in the stellar population. This is corrobo- 
rated by our study, in the sense that a much better fit is 
obtained, if one allows for a variation of the metallicity. 




age (Gyr) 

Fig. 9. Evolution of the gas surface density in the solar 
neighbourhood, computed with our SFH (solid line). The 
two dotted lines indicate the error at the ±lcr level. The 
circle at Gyr is the present surface gas density (Rana, 
1992). The circle at 15 Gyr is the present dynamical den- 
sity (Kuijken & Gilmore, 1991). For both values 1 a error 
bars are given 



The derived relation between age and metallicity is 
compared in Fig. |l^ with those obtained by others. The 
metallicity increases with time, reaching a peak (at so- 
lar metallicity) at an age of 4 Gyr; thereafter it descends 
somewhat, going through a minimum at 1 Gyr, and then 
rising again. Though Fig. [7] shows that this wavy structure 
is significant, we prefer not to place too strong an empha- 
sis on it, because on one hand the information is used 
for both SFH and AMR, that are treated as independent 
functions. With this greater amount of freedom, a good fit 
may well be achieved. On the other hand, the finding that 
the peak SFH at an age of 1.5 Gyr is not associated with 
or closely followed by a strong increase in the metallicity, 
as would be expected from chemical evolution, points well 
in the same direction. Also, the time resolution in the SFH 
is rather coarse (perhaps 0.3 dex) which also pertains to 
the AMR. 

The causal link betwen star formation and metal en- 
richment could and should be used as a constraint to make 
the SFH and AMR consistent with each other. While it 
is conceivable that it be incorporated in our method, it is 
not the aim of the present paper to do that. 
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Fig. 10. Age-metallicity relations obtained with the in- 
version (full line) in comparison with the relation of 
Meusinger et al. (1991) - dashed lines giving the aver- 
age and ±lcr curves - Crosses are individual stars from 
Edvardsson et al. (1993). 



The histogram of the metallicities of the G dwarfs is a 
powerful constraint for the chemical evolution of the disk 
(cf. Pagel, 1997). From our SFH and the time derivative 
of our AMR we compute this distribution function as 



dN*/dlgZ = 



dN^/dt 
dig Z/dt 



oc Z- i>(t(Z))/{dZ/dt) 



(13) 



which is shown in Fig. [Tl] in comparison to the histogram 
from Rocha-Pinto & Maciel (1997). Note that the distri- 
bution from Rocha-Pinto & Maciel suffers from the addi- 
tional dispersion inherent to the measurements, that we 
did not try to simulate. Our distribution has a sharp peak 
near solar metallicity, because our AMR is rather flat at 
an age around 1 Gyr. This could be due to errors in the 
metallicities as determined with Stromgren photometry. 
Using their sample of G dwarfs, Rana & Basu (1990) 
show that the intrinsic dispersion of metallicity reaches 
0.24 ± 0.10 dex at a given age. 

8. Conclusions 

The method presented here permits us to extract SFH, 
IMF slope and age-metallicity relations in a coherent 
way taking into account the selection function in appar- 
ent magnitude and the dispersion velocity of the stars. 
Indicators like the mean index and the resolution give a 
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Fig. 11. Metallicity distribution function for the G dwarfs 
(Rocha-Pinto & Maciel, 1997, dashed line) and computed 
from our SFH and AMR (full line) 



good idea of the quality of the results. It shows clearly 
that the mixed populations cannot be separated with an 
infinite accuracy and that the results depends on the ac- 
curacy of the evolutionary tracks. 

The observed Hipparcos HR diagram is well repro- 
duced by the model at the main sequence level while the 
horizontal giant branch is not correctly fit. This could be 
due to the poor knowledge of the parameters that con- 
strain the evolution of old stars. Also, our results could be 
biased by the presence of undetected binary stars or pe- 
culiar colour behaviour of A stars in the main sequence. 

The derived SFH shows enhanced star formation about 
1 to 2 Gyrs ago, while the present SFR is about equal to 
the past average SFR. This agrees well with the overall 
gas consumption in the solar neighbourhood. The AMR 
displays a systematic rise in metallicity with time. Both 
the AMR and SFH are fairly consistent with the paucity 
of metal-poor G dwarfs. 

The principal result found in this study is that the 
thin disk of the Milky Way is rather young and has a 
solar metallicity. 
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